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Abstract 

The influence of periodic and random surface textures on the flow structure and effective shp 
length in Newtonian fluids is investigated by molecular dynamics (MD) simulations. We consider 
a situation where the typical pattern size is smaller than the channel height and the local bound- 
ary conditions at wetting and nonwetting regions are characterized by finite slip lengths. In case 
of anisotropic patterns, transverse flow profiles are reported for flows over alternating stripes of 
different wettability when the shear flow direction is misaligned with respect to the stripe ori- 
entation. The angular dependence of the effective slip length obtained from MD simulations is 
in good agreement with hydrodynamic predictions provided that the stripe width is larger than 
several molecular diameters. We found that the longitudinal component of the slip velocity along 
the shear flow direction is proportional to the interfacial diffusion coefficient of fluid monomers in 
that direction at equilibrium. In case of random textures, the effective slip length and the diffusion 
coefficient of fluid monomers in the first layer near the heterogeneous surface depend sensitively 
on the total area of wetting regions. 

PACS numbers: 68.08.-p, 83.50.Rp, 47.61.-k, 83.10.Rs 
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I. INTRODUCTION 



Modeling fluid flows over chemically or topographicallypatterned substrates is important 
for micro and nanofluidic applications involving mixing |l| and separation processes As 
the surface to volume ratio increases, the role of hydrodynamic boundary conditions in 
determining fluid velocity profiles becomes dominant. It is well recognized now that the 
classical no-slip boundary condition can, in principle, be violated and the velocity profiles 
can be significantly affected by the interfacial slip [Sj . The degree of slip is usually quantified 
in terms of the Navier slip length, which is defined as a distance between the real interface 
and imaginary plane where the extrapolated tangential velocity component vanishes. The 
magnitude of the slip length for smooth nonwetting surfaces is typically on the order of 
tens of nanometers ^-[Sj; however, in special cases of nanoengineered super hydrophobic 
surfaces slip lengths in the micron range were reported 9l-|l2j. Although the Navier-Stokes 
equation with slip boundary conditions is often used to model small scale flows, the limits 
of validity of the continuum description of complex flows at nanometer scales remain not 
fully understood 

In recent years, a number of molecular dynamics (MD) studies have examined factors 
that determine the magnitude of the slip length at interfaces between crystalline surfaces 



and monatomic liquids 
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3l-|27j. One of the most important conclusions is that the degree 
of slip strongly correlates with the in-plane structure in the first fluid layer induced by the 
periodic surface potential 15||. An estimate of the low-shear-rate limit of the slip length 
can be obtained via the Green-Kubo relation between the friction coefficient at the interface 
and the time integral of the autocorrelation function of the lateral force that acts on the 



adjacent fluid from the solid wall 



18| . In general, when the surface energy is weak, the slip 



length is constant only at relatively low shear rates and it increases nonlinearly at high shear 



rates 



16j. It was recently demonstrated that the linear regime of slip holds when the slip 



velocity of the first fluid layer is smaller than the diffusion rate of fluid monomers over the 
distance between the nearest minima of the periodic surface potential at equilibrium 28 |. 
It was also found that at sufficiently high shear rates, the slip length becomes anisotropic 



for dense walls with weak surface energy; and, in particular, the slip 
the flow is oriented along the crystallographic axis of the wall lattice 



eng th increases when 



28|. 



Several studies considered the flow of a Newtonian fluid over surfaces patterned with 
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stripes o 



tions 
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different wettability using both molecular dynamics and continuum simula- 



-|32|. In the presence of flow, a heterogeneous surface with mixed boundary condi- 



tions induces spatial variations in the velocity profiles. The flow profiles averaged on length 
scales larger than the typical pattern size can be used to define the effective slip length, 
which describes the flow away from the surface. The comparative analysis between MD 
and continuum simulations demonstrated that there is an excellent agreement between the 
velocity profiles and effective slip lengths for stripe widths larger than approximately 30 
molecular diameters for flow configurations either parallel or perpendicular to the stripe 
orientation 3C|. Later studies have shown that similar conclusions hold for slip flows of 



Newtonian 33, 



34| and polymeric 35| fluids over periodically corrugated surfaces. In a 



more general situation, when the mean flow direction is not aligned with the symmetry axis 
of surface patterns, it is expected that the slip velocity will acquire a non-zero transverse 
component. In MD simulations, the transverse velocity profiles were reported in force-driven 
flows over flat surfaces with asymmetric distribution of wetting regions 32|| and in spiral 
flows inside a cylindrical channel 36|]. One of the goals of the present study is to perform a 
detailed comparative analysis of the effective slip length and velocity fields in flows around 



37- 



anisotropic textured surfaces using MD simulations and recent analytical results 

The Navier slip boundary condition for flows over arbitrary patterned surfaces can be 
formulated in tensor form, i.e., the apparent slip velocity vector is equal to the product of 



normal traction and an interfacial mobility tensor 40l-l42| . It was suggested that at the mi- 
croscopic level, the mobility tensor is related to the interfacial diffusivity per unit area 40 |. 
Recently, it was also proven that for steady noninertial flows over surfaces perturbed by arbi- 
trary periodic height and local slip fluctuations the mobility tensor is always symmetric 43 |. 
Using the theory of transport in heterogeneous media, rigorous bounds on the effective slip 
length were obtained for arbitrary two-component texture with given area fractions and 
local slip lengths 4^, |45|]. In particular, it was shown that parallel (perpendicular) stripe 
orientation with respect to the mean flow direction results in maximum (minimum) slip 



flow in a thin channel 



44| . In case when the mean flow is not parallel or perpendicular to 



37|, 



40| 



stripes, the angular dependence of the effective slip length was derived analytically 
and later verified by lattice Boltzmann simulations 46|]. One of the motivations of the cur- 
rent study is to examine the microscopic justification of the tensor formulation of the Navier 
slip condition for surfaces with anisotropic textures. 
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In this paper, the velocity fields and diffusion of fluid molecules near interfaces between 
simple fluids and surfaces patterned with anisotropic and random textures are studied by 
molecular dynamics simulations. For flows around parallel stripes of different wettability, the 
transverse and longitudinal velocity components and the effective slip length are compared 
against hydrodynamic predictions. We will show that the directional diffusion coefficient 
for fiuid molecules in contact with wall atoms correlates well with the effective slip length 
as a function of the flow direction with respect to the stripe orientation. In case of random 
wetting patterns, the effective slip length depends sensitively on the total area of wetting 
regions, in agreement with simple physical arguments. 

The rest of the paper is organized as follows. The details of molecular dynamics simu- 
lations, parameter values, and thermostatting procedure are described in the next section. 
The results for the effective slip length at surfaces with periodic and random textures and the 
numerical analysis of the interfacial diffusion of fluid molecules are presented in Section lllli 
The conclusions are given in the last section. 



II. DETAILS OF MOLECULAR DYNAMICS SIMULATION MODEL 



The schematic setup of the channel geometry and the orientation of wetting and nonwet- 
ting regions of the stationary lower wall are illustrated in Figured! The steady shear flow is 
induced by the upper wall moving with a constant velocity in the xy plane. The structure 
of three-dimensional flow flelds in the conflned fluid is determined by the local boundary 
conditions at the patterned lower wall and the orientation of the upper wall velocity with 
respect to the stripe direction. 

The pair interaction between fluid monomers {Nf = 6000) is modeled via the truncated 
Lennard- Jones (LJ) potential 



VLj{r) = 4e 



a 
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r 



(1) 



where e and a are the energy and length scales of the fluid phase, and the cutoff radius is 
rc = 2.5 (7. Wall atoms interact with fluid monomers through a modifled LJ potential with 
adjustable strength of the attractive term 



VLj{r)=Ae 



a\ 12 



(2) 
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where the parameter 5 = 1.0 for wetting regions and 5 = 0.1 for nonwetting regions and the 
rest of the parameters are the same, i.e., e„f = e, a^f = a and rc = 2.5 a. The wall atoms do 
not interact with each other. For the results presented in the next section, the parameter 
5 = 1.0 was fixed for the upper wall atoms, while the lower wall is either homogeneous {6 = 1.0 
or 6 = 0.1) or patterned with periodic stripes or random wetting regions. 

The solid walls are constructed of two layers of the face-centered cubic (fee) lattice with 
density pw = 2.3 cr"^. Each layer is composed of 576 lattice sites arranged on the (111) plane 
with [112] orientation parallel to the x direction. The nearest-neighbor distance between 
lattice sites within the (111) plane is 0.85a. The wall atoms are attached to the lattice 
sites by harmonic springs. The system dimensions in the xy plane (parallel to the confining 
walls) were kept fixed = 17.67 a and Ly = 20.41(7, and the distance between wall lattice 
planes in contact with fiuid molecules was set h = 21.54 a. The fiuid density is defined as 
a ratio of the total number of fiuid monomers to the volume accessible to the fiuid phase 
p = N f / LxLy{h — a) = 0.81 cr^^. Periodic boundary conditions for fiuid monomers and wall 
atoms were imposed along the x and y directions. The motion of the upper wall with 
a constant velocity oriented at an angle 9 with respect to the x axis was modelled by 
translating the fee lattice sites and applying periodic boundary conditions in the xy plane 
at each time step. Figure[2] shows a snapshot of the fiuid phase confined between atomically 
smooth walls. In this particular case, the upper wall velocity is oriented perpendicular to 
the stripe direction. 

The motion of wall atoms was coupled to an external heat bath by adding Langevin noise 
and friction terms to all three components of the equations of motion, e.g., in the x direction 
the equation is given by 



where the mass of a wall atom is = 100 m, the friction coefficient is r = 2.0r^^, and fi 
is a random force with zero mean and variance {fi{0)fj{t)) = 2mkBTT5{t)5ij determined 
from the fiuctuation-dissipation theorem. The temperature of the Langevin thermostat is 
T = l.le/kB, where ks is the Boltzmann constant. The wall atoms were tethered to the 



K = 2000 6 / cr"^ . It was previously shown that a sufficiently large value of the stiffness coef- 
ficient does not affect the slip length at the interface between monatomic fluids and dense 





fee lattice sites under the harmonic potential Vgp =^Kr 



with the spring stiffness coefficient 
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crystalline walls 23| . The equations of motion for wall atoms and fluid monomers were solved 
using the fifth-order Gear predictor-corrector scheme with a time step At = 0.005 r, 
where r= yjma'^ je is the LJ time scale. Thus, the oscillation time ^tx^vtlwIk, ^ 1.4 r of 
wall atoms is much larger than the integration time step. Typical values for liquid argon 
are o = 0.34nm, e/ke = 120 K and r = 2.16 x lO'^^ g |47|_ 

A common practice in non-equilibrium MD simulations is to apply the Langevin thermo- 
stat only in the direction of motion perpendicular to the plane of shear in order to maintain 



constant temperature of the fluid phase [l5|, [l6j. It is expected, however, that if the ori- 



entation of the upper wall velocity is restricted to 0° < 6 < 90° in the system shown in 
Fig.[T], then the fluid flow near the lower patterned wall will have non-zero components of 
the averaged velocity fields in all spatial dimensions, and, therefore, the Langevin friction 
term applied to fluid monomers might bias the flow profile. To test the dependence of our 
results on the thermostatting procedure, we performed simulations with a = 8.84 a and the 
upper wall velocity U = O.lcr/r oriented at ^ = 45° relative to the x direction for two cases 
where the Langevin thermostat with F = 1.0 r^^ was applied to the equations of motion for 
fiuid monomers in the direction either perpendicular to the plane of shear or parallel to 
the z direction. After averaging over thermal fiuctuations, we observed a slight difference 
in the velocity profiles near the lower wall which resulted in a discrepancy between slip 
lengths of about 0.5 cr. To eliminate the uncertainty associated with the friction term, in 
the present study, the Langevin thermostat was applied only to the equations of motion 
for wall atoms. The viscous heat in the fiuid phase was efficiently removed via interaction 
of fluid monomers with thermal wall atoms so that the fiuid temperature remained con- 
stant Tf = (1.10 ±0.01) e/ks even at the highest upper wall velocity U = 0.1 cr/r. A similar 
thermostatting procedure in sheared fluids conflned by thermal walls was implemented in 



previous MD studies |27l. |48|. 



In the present study, the simulations were performed at relatively low shear rates 
7 ^ 0.005 r^^, which required averaging of the velocity proflles over long time intervals 
(up to 10^ r) within horizontal bins of thickness Az = 0.01 a. Likewise, the fluid density 
proflles were computed within narrow bins of thickness Az= 0.01 a to resolve flne details of 
the layered structure near interfaces. The structure factor was computed in the flrst fluid 
layer according to >S'(k) = | X^^i ""^ P/ where k is a two-dimensional wave vector, 
Tj = {xj, Uj) is the position of the jth monomer, and A^^ is the total number of monomers 
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within the layer [l5|. The fluid viscosity fi = (2.15 ± 0.15) era ^ was previously found to be 
rate independent for 7 < 0.072 and insensitive to the temperature variation in the range 



of 1.1 ^ Ths/s ^ 1.35 j33|. An estimate of the maximum value of the Reynolds number 
based on U = 0.1 a/r and no-slip boundary conditions is Re = phU/fi ~ 0.77, which is 
clearly indicative of laminar flow in the channel. 

III. RESULTS 

A. Fluid density and velocity profiles for homogeneous walls 

We flrst consider steady shear flow in the channel with a homogeneous lower wall, which 
is either wetting {6 = 1.0) or nonwetting (5 = 0.1). It is well known that the presence of a 
flat crystalline surface promotes the formation of a layered structure in the adjacent fluid. 



e.g., see Ref. 49[. An example of fluid density proflles is shown in Fig.|3](a). As is evident, 
the fluid layering is most pronounced near interfaces, the amplitude of density oscillations 
gradually decays on distances of about flve molecular diameters away from the solid walls, 
and the fluid density is uniform in the middle of the channel. 

The interaction between fluid monomers and wall atoms is controlled by the strength 
of the attractive term in the LJ potential. When the parameter 6 in Eq. is reduced, 
then the well depth of the potential function decreases and the pairwise separation, where 
the potential energy reaches a minimum value, increases. For comparison, the minimum 
of the full LJ potential with 5 = 1.0 in Eq. ([2]) occurs at r = ~ 1.12 a and equals 

Vlj{1.12 a) = —e, while the modifled LJ potential with 6 = 0.1 has a much lower well depth 
Vlj{1.65 a) = — O.Ole. Therefore, it is not surprising that the amplitude of the flrst fluid 
layer near the nonwetting lower wall in Fig.|3](a) is signiflcantly reduced and its location 
is shifted away from the wall. Note also that, due to a pronounced fluid layering near 
the wetting lower wall, the amplitude of density oscillations near the upper wall is slightly 
smaller in the wetting case. Finally, we have checked that the density proflles reported in 
Fig.[3](a) for the upper wall velocity U= 0.1 a/r are the same as those computed in the 
absence of shear flow (not shown). 

The averaged fluid velocity proflles for wetting and nonwetting lower walls are presented 
in Fig.[3](b). In both cases, the velocity profiles are linear across the channel; however, the 
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slip velocity is much higher near the nonwetting lower wall. A slight downward curvature in 
the velocity profile near the nonwetting surface might be related to the fact that only fluid 
monomers with relatively large velocity component in the z direction can penetrate more 
deeply into the solid wall and, thus, the tangential velocity in that region is computed from 
biased velocity distribution. To compute the slip length, we define the location of liquid-solid 
interfaces (see vertical dashed lines in Fig. [3]) at the distance 0.5 cr away from the fee wall 
lattice planes in contact with fluid monomers. For flows over homogeneous walls, the slip 
length is determined from the linear fit to the velocity profiles excluding regions of about 2 a 
near interfaces. Depending on the wall-fluid interaction, the slip lengths are = (3.6±0.8) a 
for wetting surfaces and 6„ = (156± 10) a for nonwetting surfaces. Within the reported error 
bars, the slip lengths are independent of the shear flow orientation relative to the fee wall 
lattice. It was shown that at low shear rates 7 < 0.005 r^^ considered in the present study, 
the slip length at an interface between smooth crystalline walls and monatomic fluids is rate 
independent IgI ]. 



B. Anisotropic slip lengths for periodically patterned walls 

We next study the effects of shear flow orientation and stripe width on the flow structure 
in the channel with the lower stationary wall patterned with alternating stripes of different 
wettability, as shown in Fig.[Tl In this geometry, the textured lower wall induces wavy 
perturbations in simple shear flow which penetrate into the fluid domain on a length scale 
of about a stripe width, and the slip velocity at the lower wall is not parallel to the upper 
wall velocity when 0° < < 90°. In this section, we only consider alternating stripes of 
equal width, which are measured a = L^/ n, where integer n = 2, 4, 8, 12, and 24. Thus, 
in all cases examined, the stripe width (a/cr= 8.84, 4.42, 2.21, 1.47, 0.74) is smaller than 
the channel height h = 21.54 a; and, therefore, the longitudinal velocity profiles are expected 
to be linear across the channel except in the region of about a near the lower wall. In what 
follows, we denote the longitudinal component of the fluid velocity profile (parallel to the 
upper wall velocity) by u^\{z) and the transverse component of the velocity profile by u±{z), 
which is perpendicular to the direction of U. 

The problem of anisotropic slip flow over an array of periodic stripes of mixed wettability 
was addressed analytically assuming that the stripe width is much smaller that the fluid 
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domain 37H40| . In particular, it was shown that the angular dependence of the effective 
slip length is given by 

Ls{e) = b^cos'^e + b\\stn'^e, (4) 

where b± and 6|| are respectively slip lengths for flows perpendicular {6 = 0°) and parallel 
{9 = 90°) to the stripe orientation [37|,l39|. If the local shp lengths at wetting and nonwetting 
regions are finite and independent of the flow direction, then b\\ > b± and in the special case 
of stick-perfect slip stripes bu 



mixed boundary conditions 



2bi 



39 



38 



50| . From the solution of the Stokes equation with 



39|, the ratio of slip velocities in longitudinal and transverse 



(5) 



directions can be calculated as a function of the flow orientation 

u]_ {b\\ — b±) sin9 cos 9 
M^l bj^cos'^9 + b\\siTn?9 

In the limiting cases when the flow direction is perpendicular or parallel to the stripe ori- 
entation or when the surface is homogeneous (i.e., b\\ = b±), the transverse slip velocity 
in Eq. (|5]) vanishes. As an aside, the continuurn analysis also predicts that the transverse 
velocity component is maximum when 9 = 45° 



Examples of longitudinal and transverse velocity proflles for the smallest (a = 0.74 cr) 
and largest (a = 8.84 cr) stripe widths are presented in Fig.Hlfor selected values of 9. In 
both cases, the longitudinal velocity component is maximum (minimum) when the upper 
wall velocity is parallel (perpendicular) to the stripe orientation, and the transverse flow 
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If the 



is maximum when 9 = 45°, in agreement with the continuum analysis 
stripe width is about the molecular size, the alternating surface potential [S = 1.0 and 0.1 in 
Eq- (E]) ] represents an effectively roughened surface for the flow component perpendicular to 
the stripe orientation. As a consequence, the location of the flrst fluid layer varies periodically 
above wetting and nonwetting regions [e.g., see Fig.|3](a)], and the slip velocity along the 
X direction (perpendicular to stripes) is reduced. In contrast, when the flow is parallel to 
stripes, fluid monomers are transported along homogeneous wetting or nonwetting regions, 
and the effect of surface roughness is absent. This explains the relatively large variation 
of the longitudinal slip velocity as a function of the flow orientation for the stripe width 
a = 0.74 cr. For the largest stripe width a = 8.84 cr, the velocity proflles acquire pronounced 
oscillations near the lower wall because of the mismatch between the location of peaks in 
density proflles above wetting and nonwetting regions [shown in Fig.|3](a)]. Similar effects 
were reported in the previous study where only flows parallel and perpendicular to stripes 



were considered 



30|. For the results presented below, the effective slip length was computed 



by extrapolating the linear part of the longitudinal velocity profiles to zero velocity. 

The angular dependence of the effective slip length for the indicated values of the stripe 
width is presented in Fig.|3 As expected, Lg monotonically increases as 9 approaches 90°. 
For a given stripe width, the values 6_l and h\\ were determined from the longitudinal velocity 
profiles and then used in Eq. (jlj) to compare the results of MD simulations with continuum 
predictions (shown by red curves in Fig. [5]). The agreement between the MD data and 
continuum solution Eq. (jlj) becomes progressively better as the stripe width increases up 
to a = 8.84 cr. It should be mentioned that these results do not contradict the conclusion 
drawn from the previous study 30|], which demonstrated that the agreement between con- 
tinuum analysis and MD simulations holds when the stripe width is larger than about thirty 
molecular diameters. In the context of the present study, this conclusion could be verified 
by computing Ls{9) directly from the solution of the Stokes equation for the flow geometry 
shown in Fig.[T]with the local slip lengths at wetting and nonwetting regions extracted from 
MD simulations. Such an analysis was not performed in the current study. 

As discussed above, the transverse flow appears when the upper wall velocity is neither 
parallel nor perpendicular to the stripe orientation. Next, we present a more detailed com- 
parative analysis of the transverse slip velocity at the lower wall based on the MD data 
and the continuum solution Eq. ([5]). In MD simulations, the longitudinal and transverse slip 
velocity components of the first fluid layer were computed as follows: 

■"1.11= / 'u±^\\{z) p{z) dz I I p{z)dz, (6) 

where the limits of integration [zq = — S.lcr and zi = —5.8 a) are defined by the width of 
the first peak in the density profile above the patterned lower wall. The ratio u]_/u^^^ as 
a function of 9 is plotted in Fig.O The results show that although the magnitude of the 
transverse slip velocity is largest when 9 = 45°, the maximum angle between the upper wall 
velocity and the slip velocity occurs when 9 < 45°. For example, this angle is about 25° for 
the stripe width a = 1.47 a when 9 = 30°. Further, the MD values bj_{a) and b\\{a) were 
used in Eq. ([5]) to compute the ratio u]_/u^^^ as a function of 9 (see red curves in Fig.|6]). 
The continuum solution Eq. ([5]) is only in qualitative agreement with the MD data. The 
discrepancy might be attributed to the roughness effect discussed earlier and to the fact that 
the stripe width is comparable to the fluid molecular size. Another possible contributing 
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factor is the uncertainty in defining the exact location of the hquid-sohd interface used for 
computing the effective shp length (e.g., see the vertical dashed lines in Fig.H]). Remember 
that the position of the first fluid layer is displaced by about a from the fee lattice plane. 
Thus, if the location of the interface is taken at the position of the first fluid layer, then 
the effective slip lengths become — t- || +0.5 a. The corresponding continuum solution 
Eq. (|5]) is shown by the blue dashed curve in Fig. |6] when the stripe width is a = 4.42 a. 
Notice the small difference between continuum solutions, which, in general, is expected to 
be negligible when 6^ || ^ 0.5 cr. For completeness, the results for the smallest stripe width 
a = 0.74(7 are reported in Fig.[71 In this case, the discrepancy between the MD data and 
continuum solutions Eqs. dl])-® is most pronounced. 



C. Interfacial diffusion near surfaces of patterned wettability 

A tensorial generalization of the Navier slip condition for flows over anisotropic surfaces 
involves a relation between the normal traction at the interface and fluid slip velocity via 



an interfacial mobility tensor 40|]. In analogy with the theory of Brownian motion, it was 



conjectured that the mobility of fluid molecules near anisotropic surfaces is directly related 



to the interfacial diffusivity per unit area j40|. Simply put, it implies that the effective slip 



length at the interface between a Newtonian fluid and a textured surface is proportional to 
the diffusion coefficient of fluid molecules near the surface. In the present study, the diffusion 
coefficient was estimated from two-dimensional trajectories of fluid monomers within the first 
layer near the patterned wall at equilibrium (i.e., when both walls are at rest). 

The numerical analysis of the molecular displacement was performed only for those fluid 
monomers that remained in contact with the lower wall atoms (in the first fiuid layer) 
during the diffusion time interval. It should be noted that it is important not to include 
fluid monomers further away from the patterned wall because their diffusion in the xy plane 
quickly becomes isotropic. For example, it was recently shown that when the LJ fluid is 
confined in narrow slit-pores (with the channel height of about 5 a) and both walls are 
patterned with stripes of different wettability, then the mean square displacement curves, 
which were averaged over all fluid molecules in the direction either parallel or perpendicular 



to stripes, nearly coincide with each other 



5l|. 



FigurelH] shows the time dependence of the mean square displacement (MSD) of fluid 
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monomers near the lower wall with stripes of width a = 2.21 o". For each curve plotted in 
Fig. El the displacement vector in the xy plane was projected onto a line oriented at an angle 
^ = 0°, 25°, 45°, 60°, 90° with respect to the x axis. We find that the diffusion is isotropic 
only when the average displacement of fiuid monomers is less than the stripe width. When 
tq > a, the mean square displacement curves exhibit a gradual crossover to a linear regime 
(tq ~ t) where the diffusion becomes anisotropic. Averaging over long time periods (up 
to 2.6 X 10^ r ^ 5.6 x 10^^ s) was required to accumulate good statistics because of the 
relatively low probability that a fluid monomer would remain in the flrst layer for a long 
time; especially above the nonwetting regions where the fluid density layering is reduced and 
fluid monomers jump in and out of the flrst layer more frequently. 

The diffusion coefficient was computed from the linear slope of the mean square displace- 
ment as a function of time (r^ = ADe t) in the regime t > 145 r, which is shown in more detail 
in the lower inset of Fig. [HI The slight nonlinearity in MSD curves at large times is refiected 
in the error bars for the diffusion coefficient. Most interestingly, as shown in the upper inset 
of Fig.lHl there is an almost linear correlation between the in-plane diffusion coefficient and 
the effective slip length as a function of 6. It means that the longitudinal component of 
the slip velocity along the shear flow direction is proportional to the diffusion rate of fluid 
monomers along that direction at equilibrium. The results in Fig. [8] provide support for the 
microscopic justification of the tensorial formulation of the Navier slip boundary condition 
in the case of Newtonian fluids and molecular-scale surface textures. A similar correlation 
between Dg and Ls also holds for smaller stripe widths (a = 0.74 a and a = 1.47 cr), although 
the difference in diffusion rates along = 0° and 90° directions becomes smaller when the 
stripe width decreases (not shown). For larger stripe widths (a = 4.42 cr and a = 8.84 a), 
an accurate resolution of the mean square displacement curves in the linear regime {re > a) 
would require a long averaging time because of the low probability that a fluid monomer 
will remain within the flrst layer for a long time interval. 

D. Slip flows over surfaces with random textures 

In general, the problem of slip flow over surfaces with mixed boundary conditions specifled 
on randomly distributed regions is difficult to treat either analytically or numerically. An 
example of a statistical analysis of the effective slip boundary condition for liquid flow over 
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a plane boundary with randomly distributed free-slip regions was presented in Ref. 52 1. 



It was found that the effective slip length is proportional to the typical size of free-slip 



regions and a factor that depends on the fractional area coverage [52[. In this section, we 
consider a simple shear flow over a smooth substrate with random distribution of wetting 
and nonwetting regions. The system setup is essentially the same as described in Sec.HTl 
except that the parameter 6 in Eq. (|2]) is randomly chosen to be either 0.1 or 1.0 for the 
lower wall atoms. Each of two fee lattice layers of the lower wall contains the same number 
of weakly {6 = 0.1) or strongly {6 = 1.0) attractive atoms. In what follows, a fraction of 
wall atoms with 6= 1.0 is denoted by (p. Due to limited computational resources only one 
realization of disorder was considered for each value of 0. 

FigurelH] shows the effective slip length as a function of (p for four orientations of the upper 
wall velocity relative to the lower wall. As expected, Lg decreases with increasing the total 
area of wetting regions. Note that for each value of 0, the flow is almost isotropic; the slight 
discrepancy in most probably due to flnite size effects. When = 0.5, the averaged effective 
slip length Lg ~ 7.0 a is between b± and b\\ for any value of the stripe width reported in 
Figs.[5]and[71 which conflrms earlier conclusions that parallel (perpendicular) stripes attain 
maximum (minimum) slippage. In the limiting cases — )■ or — t- 1, the fluid velocity 
flelds near the lower wall are parallel to the xy plane, and the total friction coefficient (the 
ratio of fluid viscosity to slip length) can be estimated by simply adding contributions from 
wetting and nonwetting follows: 

// /i ^ /^(l - 0) ^^-j 



Ls{(f)) bn 

where b^ and 6„ are respectively slip lengths for wetting (0 = 1) and nonwetting (0 = 0) 
surfaces. This immediately gives the effective slip length as a function of 

As shown in Fig.[9l the agreement between the MD data and Eq. ([8]) is quite good for all 
in the range [0, 1]. However, this correspondence might, in general, not hold at intermediate 
values of and larger system size in the x and y directions because of the spatial variation 
of velocity proflles induced by the heterogeneous surface. Interestingly, the formula for 
the effective friction coefficient, Eq. ([7]), also accurately describes hydrodynamic flows along 



alternating stripes with local slip lengths that are larger than the system size 
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53]. In 



addition, it was shown numerically that an interpolation formula like Eq. ([7]) predicts the 
effective slip length for composite interfaces which consist of periodically distributed solid 
and gas areas [54 1. 

Similar to the analysis of the interfacial diffusion presented in Sec. 1111 C| we evaluate the 
in-plane diffusion coefficient in the absence of shear flow for fluid monomers in the first layer 
near the lower wall for the same realization of disorder as in Fig.|9l As an example, a typical 
trajectory projected onto the xy plane is shown in the inset of Fig.[TOl It can be seen that 
when = 1, the diffusive motion of a fluid monomer in contact with lower wall atoms is 
strongly influenced by the periodic surface potential; most of the time the monomer resides 
near the local minima of the surface potential. In this case, the surface-induced structure 
in the first fiuid layer is quantified by a distinct peak in the structure factor at the main 
reciprocal lattice vector ^(S.SS o"~^, 0) ~ 1.7, which is comparable to the height of a circular 
ridge S{2Ti/a) ^ 2.9 characteristic of the short range order. As the fraction of wetting 
regions decreases, the amplitude of density oscillations near the lower wall is reduced [see 
Fig.[3](a)], trajectories of fiuid monomers become less affected by the corrugation of the 
random surface potential, and the time interval between jumps in and out of the first fiuid 
layer decreases. When (p < 0.08, the in-plane structure factor, estimated in the first fiuid 
layer at equilibrium, does not contain any peaks at the reciprocal lattice vectors. 

The mean square displacement curves as a function of time are displayed in Fig. [TO] for 
selected values of (p. It is apparent that the diffusion becomes faster as the total area 
of nonwetting regions increases. The in-plane diffusion coefficient was estimated from the 
Einstein relation r'^y = ADxyt when t > 6t. As shown in Fig.[TTl the diffusion coefficient 
gradually varies between two values obtained for homogeneous surfaces with = and 
= 1. Furthermore, a correlation between the effective slip length (averaged over four 
orientations of the mean flow) and the in-plane diffusion coefficient is presented in the inset 
Fig.[TTl These data indicate a nearly linear dependence between Lg and D^y when the 
fraction of wetting regions is large; and, as discussed earlier, there is a strong coupling 
between the diffusion of fluid monomers in the flrst layer and the periodic surface potential. 
In the opposite limit of small 0, the net adsorption energy is reduced and the in-plane 
diffusion of fluid monomers is mostly dominated by the interaction with its fluid neighbors. 
In this regime, the effective slip length increases rapidly as the fraction of wetting regions 
decreases. 
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IV. CONCLUSIONS 



In this study, molecular dynamics simulations were performed in order to investigate 
the effective slippage and molecular diffusion at surfaces patterned with periodic or random 
textures. In our setup, the typical size of surface patterns is smaller than the channel dimen- 
sions, and the local boundary conditions at homogeneous wetting or nonwetting surfaces are 
described by finite slip lengths. Particular attention was paid to the implementation of a 
thermostatting procedure that does not bias flow profiles and diffusion of fluid monomers. 

For flows over surfaces patterned with stripes of different wettabihty, the heterogeneous 
surfaces induce wavy perturbations in velocity profiles and the slip velocity acquires a trans- 
verse component. In this case, the effective slip length depends of the shear flow direction 
with respect to the stripe orientation. We found that the angular dependence of the effec- 
tive slip length computed by molecular dynamics simulations agrees well with the analytical 
solution of the Stokes equation provided that the stripe width is larger than several molec- 
ular diameters. At the same time, however, the ratio of the transverse and longitudinal 
components of the slip velocity agrees only qualitatively with hydrodynamic predictions. 
Furthermore, the interfacial diffusion coefficient of fluid molecules correlates well with the 
effective shp length as a function of the shear flow direction. The numerical analysis was 
performed only for fluid monomers that remain in contact with the wall atoms during the 
diffusion time interval. These findings lend support for the microscopic justification of the 
tensor formulation of the effective slip boundary conditions for noninertial flows of Newto- 
nian fluids over smooth surfaces with nanoscale anisotropic textures. 

In case of random surface textures, the simulation results and simple physical arguments 
show that the effective slip length is determined by the total area of wetting regions. When 
the fraction of wetting regions is large, the diffusive motion of ffuid monomers is strongly in- 
ffuenced by the periodic surface potential; and the effective slip length is nearly proportional 
to the in-plane diffusion coefficient at equilibrium. In the opposite limit of small wetting 
areas, the diffusion of fluid monomers is less affected by the corrugation of the surface po- 
tential, and the effective slip length depends sensitively on the number of strongly attractive 
wall atoms. 
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FIG. 1: (Color online) A schematic representation of the channel geometry and surface patterns 
indicated by the blue color (wetting regions) and white color (nonwetting regions). Steady shear 
flow is induced by the upper wall moving with a constant velocity U in the xy plane at an angle 9 
with respect to the x direction. The lower patterned wall is stationary. 
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FIG. 2: (Color online) A snapshot of 6000 fluid monomers (open blue circles) confined between 
atomistic walls patterned with wetting (filled gray circles) and nonwetting (filled yellow circles) 
regions. The width of stripes at the lower stationary wall is a = 4.42 a and the upper wall velocity 
is [/ = O.lcr/r. 
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FIG. 3: (Color online) Averaged (a) density and (b) velocity profiles across the channel for wetting 
6=1 (red circles) and nonwetting 5 = 0.1 (black squares) lower walls. The upper wall velocity in the 
X direction is U = 0.1 a /t. The vertical axes at z/a = — 6.67 and 14.87 coincide with the location 
of the fee lattice planes in contact with fluid molecules. The vertical dashed lines at z/a = — 6.17 
and 14.37 indicate the location of liquid-solid interfaces. 
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FIG. 4: (Color online) Averaged longitudinal (upper set of curves) and transverse (lower set of 
curves) velocity profiles for the indicated values of 9 and stripe widths a = 0.74 o" (left panel) and 
a = 8.84 cr (right panel). The vertical axes coincide with the location of the fee lattice planes at 
z/a= — 6.67 and 14.87. The vertical dashed lines at z/a= — 6.17 and 14.37 indicate liquid-solid 
interfaces. 
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FIG. 5: (Color online) Variation of the effective slip length Lg/cr as a function of 9 for the indicated 
values of the stipe width. The MD data are shown by open circles. The red curves are hydrodynamic 
predictions computed using Eq. (jl]). 
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FIG. 6: (Color online) The ratio of transverse and longitudinal components of the slip velocity 
as a function of 9. The MD data are shown by open circles. The solid red curves are continuum 
predictions calculated using Eq. ([5]). The dashed blue curve is Eq. ([5]) where y are defined with 
respect to the location of the first fluid layer (see text for details). 
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FIG. 7: (Color online) The effective slip length Lg/a as a function of 6 for the stripe width 
a = 0.74(7. The continuum prediction, Eq. is shown by the red curve. The inset shows the 
angular dependence of the ratio of slip velocities in transverse and longitudinal directions. The red 
curve in the inset is the prediction of Eq. ([5]) . 
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FIG. 8: (Color online) The mean square displacement of fluid monomers in the first layer near the 
lower patterned wall with the stripe width a = 2.21 a, when both walls are at rest. Each curve is 
computed by taking the component of the displacement vector along a line oriented at an angle 
9 = 0°, 25°, 45°, 60°, 90° with respect to the x axis (from bottom to top). The dashed line 
with unit slope is plotted for reference. The lower inset shows an expanded view of the same data 
at large t. The upper inset shows a correlation between the in-plane diffusion coefficient and the 
effective slip length for the same set of 9. The straight blue line is the best fit to the data. 
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FIG. 9: (Color online) The effective slip length L^ja versus the fraction of randomly distributed 
attractive atoms {b = 1.0) at the lower wall. The orientation of the upper wall velocity is specified 
in the inset. The black dashed curve is Eq. ([8]) with 6^(0 = 1) = 3.6 a and 6n(0 = 0) = 156 cr. 
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FIG. 10: (Color online) The mean square displacement of fluid monomers in the first layer near 
the lower wall with = 0, 0.08, 0.58, 0.75, 1.0 (from left to right) when [/ = 0. The dashed line 
with unit slope is plotted as a reference. The inset shows a typical trajectory of a fluid monomer 
for about 100 r near the wetting wall cj) = 1.0. The positions of the fee lattice atoms in the xy 
plane are denoted by open circles. 
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FIG. 11: (Color online) The in-plane diffusion coefficient Dxy (in units cP' jr) as a function of 
The straight blue line is the best fit to the data for > 0.2. The inset shows a correlation between 
the effective slip length and the diffusion coefficient for ^ (/> ^ 1. The solid curve is a guide for 
the eye. The dashed line with unit slope is drawn for reference. 
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